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INTRODUCTION 

Two-dimensional  (2D)  material  heterostructures  offer  novel  and  compelling  electronic  and 
optical  properties.  Density  functional  theory  (DFT)  is  often  used  to  derive  the  electronic  band 
structure  of  the  2D  heterostructures  from  first  principles  as  well  as  to  validate  experimental  results. 
However,  the  implementation  of  DFT  requires  an  in-depth  understanding  of  the  geometric  properties 
of  the  system  being  analyzed.  The  creation  of  a  unit  cell  that  accurately  describes  the  system 
remains  one  of  the  largest  challenges  for  DFT  calculations.  As  the  unit  cell  size  increases, 
computational  requirements  increase  exponentially.  However,  too  small  of  a  unit  cell  may  fail  to 
represent  short-  to  mid-range  order  in  the  lattice.  Balancing  these  conflicting  requirements  of  cost 
and  accuracy  is  still  one  of  the  biggest  barriers  to  implementing  DFT  calculations.  An  algorithm  to 
automatically  optimize  the  process  would  greatly  streamline  the  method  of  creating  a  unit  cell.  While 
many  methods  have  undoubtedly  been  created  for  matching  lattice  constants  of  dissimilar 
nanomaterials,  very  few  are  actually  covered  explicitly  in  literature.  To  rectify  this,  and  to  ensure 
other  researchers  will  not  have  to  resolve  this  problem,  this  report  will  review  the  underlying 
geometry  of  popular  2D  materials,  draw  comparisons  to  other  common  unit  cells,  then  offer  a 
solution  which  facilitates  the  creation  of  an  optimal  unit  cell  for  any  heterogeneous  structure 
consisting  of  multiple  layers  of  2D  and  close-packed  materials. 


TWO-DIMENSIONAL  MATERIAL  GEOMETRY  AND  ANALOGS  WITH  CLOSE-PACKED  SYSTEMS 

The  vast  majority  of  2D  materials  currently  under  research  focus,  including  graphene, 
graphane,  silicine,  germanane,  and  all  transition  metal  dichalcogenides  (TMDCs),  form  hexagonal 
planes  (ref.  1).  These  structures  are  all  easily  expressed  using  60  deg  monoclinic  Bravais  lattices 
where  the  rhombic  sides  of  the  Bravais  lattices  fit  together  to  form  hexagonal  monolayers.  When 
metals  are  deposited  as  electrodes  on  top  of  these  2D  layers,  the  face-centered  cubic  (fee)  or 
hexagonal  close-packed  (hep)  lattice  structure  of  the  metals  can  also  be  modeled  as  rhombic  unit 
cells  for  consistency  with  the  underlying  stacked  heterostructures.  The  hep  structures  are  described 
using  the  primitive  hexagonal  lattice,  with  a=  b*  c,  a  =  (B  =  90°,  y  =  60°  and  two  atoms,  at  (0,0,0) 
and  (%,%,1/2).  We  can  also  represent  a  fee  crystal  using  a  rhombic  unit  cell  by  rotating  the  unit  cell 
45  deg  and  making  the  close-packed  plane  the  base  of  the  unit  cell,  with  a  -  b±  c,  a  =  (B  =  90°,  y  = 
60°  and  four  atoms,  at  (0,0,0),  (14,14,0),  (14,0, 14),  and  (0, 14,14).  As  the  primary  difference  between 
hep  and  fee  structures  lies  in  differing  stacking  sequences  (ABAB  and  ABCABC,  respectively),  it 
makes  sense  that  a  fee  cell  taken  using  the  close-packed  plane  as  the  base  is  similar  to  the  hep  unit 
cell  but  with  one  more  layer  added. 


MATCHING  SYSTEM  LATTICE  VECTORS:  AN  OPTIMIZATION  PROBLEM 

Now  that  we  are  able  to  express  various  crystal  structures  as  hexagonal  lattices  with  a 
rhombus  base,  we  can  address  the  core  issue  of  matching  differing  lattices  into  a  larger  primitive 
cell.  In  order  for  each  of  the  materials  in  the  system  to  repeat  correctly,  the  overall  system  cell  must 
be  represented  as  a  linear  combination  of  each  individual  material’s  unit  cell  (fig.  1).  We  define  a 
supercell  for  each  material  with  lattice  constant 

lAscI;  =  Jaf  *  (mf  +  nf  -  mt*  n *)  (1) 

where  a*  is  the  lattice  parameter  for  a  given  material,  and  m;  and  n*  are  integers  representing  the 
number  of  unit  cells  in  each  direction.  We  then  seek  the  minimum  value  |asc|j  such  that  it  is  equal  for 
each  material.  By  generating  a  value  for  |asc|j  for  each  material  in  the  system  then  minimizing  the 
variance  and  magnitude,  an  amorphous  geometric  problem  turns  into  a  fairly  simple  optimization 
problem. 
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Note:  In  this  example,  m1  =m2  =  4,  nt  =  2 ,n2  =  1.  a1/a2  ~  1.04,  asc/a2  *  3.6.  Purple  spheres  denote  Mo,  Green  Se, 

and  Yellow  S. 

Figure  1 

Matching  the  system  lattice  vector  for  two  different  transition  metal  dichalcogenides  (TMDs),  MoSe2 

and  MoS2 

In  order  to  solve  the  optimization  problem,  constraints  are  necessary,  and  nt  must  be 
integers,  and  the  unit  cell  lattice  parameters  of  each  material  were  allowed  to  vary  by  up  to  1%, 
which  significantly  decreased  the  cell  size  while  imposing  minimal  total  strain  in  the  structure.  A 
maximum  value  of  6  was  set  for  rrii  and  nh  as  a  larger  system  would  prove  extremely 
computationally  expensive.  A  minimum  of  1  for  and  0  for  forces  the  system  lattice  parameter 
to  have  a  non-zero  value  but  allows  for  the  smallest  possible  system.  For  this  work,  a  genetic 
algorithm  was  used  to  solve  the  optimization  problem.  Genetic  algorithms  begin  with  a  large 
population  of  potential  solutions  and  then  undergo  an  iterative  process  mimicking  evolution.  Each 
solution  is  evaluated  compared  to  a  ‘fitness  function’  that  defines  the  end  goal  of  the  optimization 
problem.  For  this  problem,  the  fitness  function  minimized  the  variance  in  the  vector  of  each 
individual  material  asc.  Genetic  algorithms  converge  quickly  and,  depending  on  the  initial  population, 
effectively  find  global  minima.  Reference  2  provides  an  excellent  primer  on  genetic  algorithms. 

Komsa  and  Krasheninnikov  (ref.  3)  took  a  similar  approach  for  optimizing  the  unit  cell  for  a 
heterostructure  of  two-dimensional  metal  chalcogenides  (TDMCs).  However,  the  method  used  in 
this  report  differs  from  theirs  in  a  few  key  respects.  For  example,  in  considering  the  bilayer 
heterostructure  of  TDMCs  with  one  member  of  the  bilayer  being  always  MoS2  and  the  other  varied 
between  WS2,  MoSe2,  MoTe2,  BN,  or  graphene,  Komsa  and  Krasheninnikov  fixed  the  lattice 
parameter  of  MoS2  and  only  allowed  the  second  layer  to  strain.  We  let  both  of  the  bilayers  undergo 
a  strain  up  to  1  %,  allowing  us  to  more  effectively  find  a  minimum  and  to  distribute  the  strain  from  the 
lattice  mismatch  among  both  layers  of  the  heterostructure.  Additionally,  while  Komsa  and 
Krasheninnikov  chose  the  smallest  unit  cell  with  <1%  strain,  we  used  a  weighting  algorithm  to 
choose  a  unit  cell  with  near-optimal  size  and  minimal  strain. 
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GENERATING  THE  SYSTEM  UNIT  CELL 


The  found  values  for  asc,  ait  mt,  and  allow  us  to  generate  the  crystal  lattice  for  each 
material.  The  angle  a  between  the  system  unit  cell  lattice  vector  and  the  material  lattice  vector  is 
given  by 


«i  = 


cos  1(a|c  +  aj 


2  CLct-TTli  CLi 


)■ 


(2) 


The  periodicity  of  the  material  lattices  makes  generating  the  structure  unit  cell  fairly  straightforward. 
From  the  top  left  corner  of  the  system  unit  cell,  every  atom  in  a  given  material  repeats  at  every  linear 
combination  of  moving  at  in  the  a  direction  and  moving  at  in  the  a  +  60°  direction  (fig.  2). 


Figure  2 

Matching  the  system  lattice  vector  to  our  two  respective  crystal  lattices 


TRANSITION  METAL  DICHALCOGENIDES  (TMDCS)  WITH  MISMATCHED  LATTICE 

PARAMETERS 

A  bilayer  of  MoS2  and  MoSe2was  chosen  as  the  first  system  to  evaluate.  In  the  calculations, 
a  relaxed  MoSe2  layer  has  a  lattice  parameter  of  3.28  A  and  a  relaxed  MoS2  layer  has  a  lattice 
parameter  of  3.16  A.  This  genetic  algorithm  successfully  matched  the  two  layers,  as  shown  in  figure 
3.  We  also  used  it  to  match  a  series  of  different  TMDCs.  The  results  are  summarized  in  table  1 .  The 
strain  (e)  in  the  individual  layers  as  shown  in  table  1  is  less  than  1%  for  all  the  cases  considered. 
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Figure  3 

The  generated  unit  cell  of  a  MoS2/MoSe2  heterostructure  bilayer 

Table  1 

Multiple  matched  monolayers 


m1 

ni 

m2 

n2 

Strained  [A] 

(Un)Strained  a2[A] 

«[°] 

Q'SC  [A] 

S] 

£2 

MoS2/MoSe2 

4 

1 

4 

2 

3.16 

(3.28)  3.29 

16.1 

11.4 

0.00% 

0.30% 

MoS2/MoTe2 

3 

3 

3 

1 

3.14 

(3.55)  3.56 

40.9 

9.42 

0.63% 

0.28% 

M0S2/WS2 

1 

0 

1 

0 

3.18 

(3.21)  3.18 

0 

3.18 

0.63% 

0.93% 

MoS2/WSe2 

4 

1 

4 

2 

3.18 

(3.33)  3.31 

16.1 

11.48 

0.63% 

0.60% 

MoSe2/MoTe2 

5 

3 

4 

4 

3.28 

(3.55)  3.57 

23.4 

11.46 

0.00% 

0.56% 

MoSe2/WS2 

4 

2 

3 

4 

3.31 

(3.21)  3.18 

13.9 

11.46 

0.90% 

0.93% 

MoSe2/WSe2 

1 

0 

1 

0 

3.31 

(3.33)  3.31 

0 

3.3 

0.90% 

0.60% 

With  the  optimized  unit  cell  in  place,  DFT  calculations  to  calculate  the  band  structure  and 
projected  density  of  states  for  the  material  can  begin.  This  simulation  was  run  in  SIESTA  using  the 
Perdew-Burke-Ernzerhof  (PBE)  functional  for  the  exchange  correlation  contributions.  All  ions  were 
relaxed  to  less  than  0.01  eV/A.  The  initial  interlayer  spacing  was  set  to  12.5  A,  and  then  the 
structure  was  relaxed  until  reaching  an  energy  minimum  at  a  spacing  of  6.4  A.  An  8x8x1  Monkhorst- 
Pack  k-point  mesh.  The  basis  set  was  double  zeta  polarized,  and  a  mesh  cutoff  of  300  Rydberg  was 
implemented. 
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DENSITY  FUNCTIONAL  THEORY  RESULTS 

As  figure  4  shows,  the  direct  bandgap  at  the  K  point  was  calculated  to  be  1 .49  eV.  Previous 
DFT  calculations  (refs.  3  and  4)  reported  a  value  of  ~1 .25  eV  for  the  direct  bandgap  at  the  K  point  for 
the  heterostructure.  The  difference  may  be  attributed  to  the  fact  that  we  have  let  both  layers  of  the 
heterostructure  to  strain  while  Komsa  and  Krasheninnikov  (ref.  3)  allowed  only  one  layer  to  stretch  or 
compress.  Interestingly,  the  direct  bandgap  of  1 .49  eV  in  our  studies  is  comparable  to  that  of  the 
individual  layers  (ref.  4).  Consistent  with  the  previous  DFT  calculations,  it  can  be  observed  in  figure 
4  that  the  valance  band  maximum  occurs  at  the  T  point  of  the  Brillouin  zone  thereby  rendering  the 
minimum  energy  gap  to  be  indirect.  The  observed  value  of  1 .15  eV  for  the  indirect  gap  is 
comparable  to  that  reported  by  Komsa  and  Krasheninnikov  (ref.  3). 


Band  Structure  of  a  WS-,  /  MoS2  Bilayer  Projected  Density  of  States 


Figure  4 

DFT  calculation  results  of  a  WS2  and  MoS2  heterostructure 

Figure  4  also  shows  the  projected  density  of  states  for  the  stacked  heterostructure.  The 
projected  density  of  states  shows  that  for  MoS2  /  WS2  most  of  the  contribution  to  the  conduction  band 
comes  from  the  molybdenum  atoms.  While  the  valence  band  contribution  comes  slightly  more  from 
the  tungsten  atoms  than  the  molybdenum,  the  difference  in  contributions  in  the  projected  density  of 
states  is  too  small  to  make  a  meaningful  conclusion  about  the  exact  nature  of  the  valence  band. 
However,  the  result  suggests  the  interesting  consequence  of  a  type  II  band  alignment  of  the 
heterostructure  with  the  valence  band  maximum  and  conduction  band  minimum  occurring  in  WS2 
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and  MoS2  layers,  respectively.  Thus,  in  an  excitonic  optical  transition,  the  electron-hole  pair  will  be 
spatially  separated.  This  will  lead  to  interesting  optical  behavior  of  the  heterostructure. 


CONCLUSIONS 

We  have  developed  a  genetic  algorithm  that  allows  us  to  rapidly  optimize  the  unit  cells  for  the 
calculation  of  the  electronic  band  structures  of  heterostructures  of  two-dimensional  (2D) 
nanomaterials.  We  applied  the  method  to  calculate  the  band  structure  of  2D  MoS2/WS2 
heterostructure,  and  our  results  are  in  good  agreement  with  that  of  the  earlier  work  of  Komsa  and 
Krasheninnikov.  The  flexibility  of  our  method  renders  it  applicable  to  any  multilayer  heterostructure 
where  the  individual  layers  can  be  expressed  using  rhombic  unit  cells.  There  are  plans  to  use  this 
technique  to  simulate  combinations  of  different  transition  metal  dichalcogenides  (TMDCs)  so  as  to 
better  understand  how  stacking  changes  the  electronic  and  optical  properties.  Further,  this  method 
will  be  used  to  model  the  electron  transport  in  a  heterostructure  that  uses  other  2D  materials  as 
metal  electrodes. 
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